R version 3.2.4 (2016-03-10) -- "Very Secure Dishes"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[R.app GUI 1.67 (7152) x86_64-apple-darwin13.4.0]

[History restored from /Users/gholliba/.Rapp.history]

> ### Jonathan D. Klingler, Gary E. Hollibaugh, Jr., and Adam J. Ramey
> ### "Don't Know What You Got: A Bayesian Hierarchical Model of Neuroticism and Ideological Uncertainty
> ### Political Science Research and Methods
> ###
> ###
> ### Main Analyses (All Figures except Figure 5; All Tables except B-1, B-2, and B-3)
> ### Description: This file creates (nearly) all tables and figures in the main paper.
> ###
> ### Note 1: Run the Preprocessing_and_Summary.R file beforehand to generate the CCES2014.Rda file.
> ### Note 2: Run the SingleCoreX.R files beforehand to generate the hierarchical model estimates.
> ### Note 3: Also run AMBoot.R beforehand to generate the bootstrapped Aldrich-McKelvey estimates.
> ### Note 4: Make sure to change the directory to the one containing this file.
>
>
> rm(list=ls())
> set.seed(1) # replicability
> 
> 
> library(grDevices)
> library(reshape2)
> library(ggplot2)
> library(RColorBrewer)
> library(rjags)
Loading required package: coda
Linked to JAGS 3.3.0
Loaded modules: basemod,bugs
> library(plyr)
> library(coda)
> library(ggmcmc)
Loading required package: dplyr

Attaching package: ‘dplyr’

The following objects are masked from ‘package:plyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

Loading required package: tidyr
> library(AER)
Loading required package: car
Loading required package: lmtest
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading required package: sandwich
Loading required package: survival
> library(basicspace)
Loading required package: tools

## BASIC SPACE SCALING PACKAGE
## 2009 - 2016
## Keith Poole, Howard Rosenthal, Jeffrey Lewis, James Lo, and Royce Carroll
## Support provided by the U.S. National Science Foundation
## NSF Grant SES-0611974

> library(texreg)
Version:  1.35
Date:     2015-04-25
Author:   Philip Leifeld (University of Konstanz)

Please cite the JSS article in your publications -- see citation("texreg").

Attaching package: ‘texreg’

The following object is masked from ‘package:tidyr’:

    extract

> library(grid)
> library(Cairo)
> library(xtable)
> 
> # load CCES data --- make sure in the correct directory first!
> load("CCES2014.Rda")
> cces_all <- cces
> cces <- cces[!is.na(cces$self_emoti) & !is.na(cces$self_consc) & !is.na(cces$self_agree) 
+              & !is.na(cces$self_extra) & !is.na(cces$self_openn) & !is.na(cces$CC421a) &
+                !is.na(cces$newsint),]
> 
> 
> # candidate names
> candNames <- c("Obama","Cruz","Clinton","Paul","Bush","Democrats",
+                "Republicans","Tea Party","Supreme Court")
> 
> 
> # not sure and skipped are observationally (and arguably substantively)
> # equivalent to each other
> cces$totalNA <- as.numeric(cces$CC334C == "Not sure") + as.numeric(cces$CC334C == "Skipped") +  # obama
+                 as.numeric(cces$CC334D == "Not sure") + as.numeric(cces$CC334D == "Skipped") +  # clinton
+                 as.numeric(cces$CC334E == "Not sure") + as.numeric(cces$CC334E == "Skipped") +  # cruz
+                 as.numeric(cces$CC334F == "Not sure") + as.numeric(cces$CC334F == "Skipped") +  # rand paul
+                 as.numeric(cces$CC334G == "Not sure") + as.numeric(cces$CC334G == "Skipped") +  # jeb bush
+                 as.numeric(cces$CC334K == "Not sure") + as.numeric(cces$CC334K == "Skipped") +  # democratic party
+                 as.numeric(cces$CC334L == "Not sure") + as.numeric(cces$CC334L == "Skipped") +  # republican party
+                 as.numeric(cces$CC334M == "Not sure") + as.numeric(cces$CC334M == "Skipped") +  # tea party
+                 as.numeric(cces$CC334W == "Not sure") + as.numeric(cces$CC334W == "Skipped")    # scotus
> 
> # need to account for those who weren't asked about Cruz/Paul
> cces$totalAsked <- 9 - (as.numeric(cces$CC334C == "Not Asked") + # obama
+                         as.numeric(cces$CC334D == "Not Asked") + # clinton
+                         as.numeric(cces$CC334E == "Not Asked") + # cruz
+                         as.numeric(cces$CC334F == "Not Asked") + # rand paul
+                         as.numeric(cces$CC334G == "Not Asked") + # jeb bush
+                         as.numeric(cces$CC334K == "Not Asked") + # democratic party
+                         as.numeric(cces$CC334L == "Not Asked") + # republican party
+                         as.numeric(cces$CC334M == "Not Asked") + # tea party
+                         as.numeric(cces$CC334W == "Not Asked")) # scotus
> 
> 
> # rescaling the B5 to be between 0 and 1
> cces$self_openn <- (cces$self_openn - 1)/6
> cces$self_consc <- (cces$self_consc - 1)/6
> cces$self_extra <- (cces$self_extra - 1)/6
> cces$self_agree <- (cces$self_agree - 1)/6
> cces$self_emoti <- (cces$self_emoti - 1)/6
> 
> 
> # making the 150k or more level consistent
> cces$faminc_recode <- cces$faminc
> cces$faminc_recode[which(cces$faminc_recode == "$150,000 - $199,999")] <- levels(cces$faminc_recode)[17] 
> cces$faminc_recode[which(cces$faminc_recode == "$200,000 - $249,999")] <- levels(cces$faminc_recode)[17] 
> cces$faminc_recode[which(cces$faminc_recode == "$250,000 - $349,999")] <- levels(cces$faminc_recode)[17] 
> cces$faminc_recode[which(cces$faminc_recode == "$350,000 - $499,999")] <- levels(cces$faminc_recode)[17] 
> cces$faminc_recode[which(cces$faminc_recode == "$500,000 or more")] <- levels(cces$faminc_recode)[17] 
> cces$faminc_recode[which(cces$faminc_recode == "$250,000 or more ")] <- levels(cces$faminc_recode)[17] 
> cces$faminc_recode <- factor(cces$faminc_recode)
> 
> 
> # total number of NAs                  
> total.na.mod.1 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
+                                                            self_agree + self_extra + self_openn, 
+                                                            data = cces, family=binomial(link="probit"))
>                                                            
> total.na.mod.2 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
+                                                            self_agree + self_extra + self_openn + 
+                                                            I(gender == "Female") + I(2014 - birthyr) + 
+                                                            I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                                            I(race == "Hispanic") + 
+                                                            I(race == "Asian" | race == "Native American" | 
+                                                              race == "Mixed" | race == "Other" | race == "Middle Eastern"), 
+                                                            data = cces, family=binomial(link="probit"))
>                                                            
> total.na.mod.3 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
+                                                            self_agree + self_extra + self_openn + 
+                                                            I(gender == "Female") + I(2014 - birthyr) + 
+                                                            I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                                            I(race == "Hispanic") + 
+                                                            I(race == "Asian" | race == "Native American" | 
+                                                              race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                                            as.numeric(educ), data = cces, family=binomial(link="probit"))
>                                                            
> total.na.mod.4 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
+                                                            self_agree + self_extra + self_openn + 
+                                                            I(gender == "Female") + I(2014 - birthyr) + 
+                                                            I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                                            I(race == "Hispanic") + 
+                                                            I(race == "Asian" | race == "Native American" | 
+                                                              race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                                            as.numeric(educ) + I(newsint == "Most of the time") + 
+                                                            I(newsint == "Don't know"), 
+                                                            data = cces, family=binomial(link="probit"))
>                                                           
> total.na.mod.5 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
+                                                            self_agree + self_extra + self_openn + 
+                                                            I(gender == "Female") + I(2014 - birthyr) + 
+                                                            I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                                            I(race == "Hispanic") + 
+                                                            I(race == "Asian" | race == "Native American" | 
+                                                              race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                                            as.numeric(educ) + I(newsint == "Most of the time") + 
+                                                            I(newsint == "Don't know") + as.numeric(faminc_recode) + 
+                                                            I(faminc_recode == "Prefer not to say"), 
+                                                            data = cces, family=binomial(link="probit"))
>                                                            
> total.na.mod.6 <- glm(cbind(totalNA, totalAsked-totalNA) ~ I(1 - self_emoti) + self_consc + 
+                                                            self_agree + self_extra + self_openn + 
+                                                            I(gender == "Female") + I(2014 - birthyr) + 
+                                                            I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                                            I(race == "Hispanic") + 
+                                                            I(race == "Asian" | race == "Native American" | 
+                                                              race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                                            as.numeric(educ) + I(newsint == "Most of the time") + 
+                                                            I(newsint == "Don't know") + as.numeric(faminc_recode) + 
+                                                            I(faminc_recode == "Prefer not to say") + 
+                                                            I(employ == "Full-time") + I(employ == "Part-time") + 
+                                                            I(employ == "Unemployed") + I(employ == "Retired"), 
+                                                            data = cces, family=binomial(link="probit"))
> 
> # coercing the models into a list for texreg
> binom.mod.list <- list(total.na.mod.1, 
+                        total.na.mod.2,
+                        total.na.mod.3,
+                        total.na.mod.4,
+                        total.na.mod.5,
+                        total.na.mod.6)
>     
> ## table 1
> # texreg for latex (commented out since it's already in the paper)                        
> # texreg(binom.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
>                           # custom.coef.names = c("Constant",
>                                                 # "Neuroticism",
>                                                 # "Conscientiousness",
>                                                 # "Agreeableness",
>                                                 # "Extraversion",
>                                                 # "Openness",
>                                                 # "Female",
>                                                 # "Age",
>                                                 # "Age^2/100",
>                                                 # "Black",
>                                                 # "Hispanic",
>                                                 # "Other Race",
>                                                 # "Education (1 = No HS; 6 = Postgrad)",
>                                                 # "High News Interest",
>                                                 # "Unknown News Interest",
>                                                 # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
>                                                 # "Income Refused",
>                                                 # "Employed Full-Time",
>                                                 # "Employed Part-Time",
>                                                 # "Unemployed",
>                                                 # "Retired"),
>                           # reorder.coef = c(2:21,1),
>                           # include.aic = FALSE, 
>                           # include.deviance = FALSE)
> 
> # screenreg for review
> # screenreg(binom.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
>                           # custom.coef.names = c("Constant",
>                                                 # "Neuroticism",
>                                                 # "Conscientiousness",
>                                                 # "Agreeableness",
>                                                 # "Extraversion",
>                                                 # "Openness",
>                                                 # "Female",
>                                                 # "Age",
>                                                 # "Age^2/100",
>                                                 # "Black",
>                                                 # "Hispanic",
>                                                 # "Other Race",
>                                                 # "Education (1 = No HS; 6 = Postgrad)",
>                                                 # "High News Interest",
>                                                 # "Unknown News Interest",
>                                                 # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
>                                                 # "Income Refused",
>                                                 # "Employed Full-Time",
>                                                 # "Employed Part-Time",
>                                                 # "Unemployed",
>                                                 # "Retired"),
>                           # reorder.coef = c(2:21,1),
>                           # include.aic = FALSE, 
>                           # include.deviance = FALSE)
> 
> # htmlreg for file
> htmlreg(binom.mod.list, file = "table1.html", digits = 3, stars = c(0.01, 0.05, 0.1),
+                           custom.coef.names = c("Constant",
+                                                 "Neuroticism",
+                                                 "Conscientiousness",
+                                                 "Agreeableness",
+                                                 "Extraversion",
+                                                 "Openness",
+                                                 "Female",
+                                                 "Age",
+                                                 "Age^2/100",
+                                                 "Black",
+                                                 "Hispanic",
+                                                 "Other Race",
+                                                 "Education (1 = No HS; 6 = Postgrad)",
+                                                 "High News Interest",
+                                                 "Unknown News Interest",
+                                                 "Income (1 = <10k; 12 = >150k; 13 = Refused)",
+                                                 "Income Refused",
+                                                 "Employed Full-Time",
+                                                 "Employed Part-Time",
+                                                 "Unemployed",
+                                                 "Retired"),
+                           reorder.coef = c(2:21,1),
+                           include.aic = FALSE, 
+                           include.deviance = FALSE)
The table was written to the file 'table1.html'.

> 
> 
> # calculating the model with the smallest BIC for plotting purposes
> binom.min.bic <- which(c(BIC(total.na.mod.1), BIC(total.na.mod.2), 
+                          BIC(total.na.mod.3), BIC(total.na.mod.4), 
+                          BIC(total.na.mod.5), BIC(total.na.mod.6)) == 
+                        min(c(BIC(total.na.mod.1), BIC(total.na.mod.2), 
+                              BIC(total.na.mod.3), BIC(total.na.mod.4), 
+                              BIC(total.na.mod.5), BIC(total.na.mod.6))))
> 
> # creating a mode function for the predicted probabilities
> Mode <- function(x) {
+   ux <- unique(x)
+   ux[which.max(tabulate(match(x, ux)))]
+ }
> 
> ## predicted probabilities for plotting (plot to be called later in combined plot)
> # "some of the time" is the median value, and a plurality of people do NOT choose "most of the time"
> # setting up the covariate values for prediction
> total.na.preds <- predict(binom.mod.list[[binom.min.bic]], 
+                           newdata = data.frame(self_emoti = seq(mean(cces$self_emoti) - 2*sd(cces$self_emoti), 
+                                                                 mean(cces$self_emoti) + 2*sd(cces$self_emoti), 
+                                                                 length.out = 100),
+                                                self_consc = mean(cces$self_consc),
+                                                self_agree = mean(cces$self_agree),
+                                                self_extra = mean(cces$self_extra),
+                                                self_openn = mean(cces$self_openn),
+                                                newsint = levels(cces$newsint)[median(as.numeric(cces$newsint))],
+                                                educ = Mode(cces$educ),
+                                                employ = Mode(cces$employ),
+                                                faminc_recode = Mode(cces$faminc_recode),
+                                                gender = Mode(cces$gender),
+                                                birthyr = median(cces$birthyr),
+                                                race = Mode(cces$race)), se = TRUE)
> 
> # predicted probabilities with 90% CIs
> # reversed scale because moving from emotional stability to neuroticism
> total.na.predvals <- data.frame(neuroticism = seq(2, -2, length.out = 100),
+                                 low = pnorm(total.na.preds$fit - qnorm(0.95)*total.na.preds$se),
+                                 mid = pnorm(total.na.preds$fit),
+                                 high = pnorm(total.na.preds$fit + qnorm(0.95)*total.na.preds$se))
> 
> 
> # standard aldrich-mckelvey
> # need to remove those where a question wasn't asked (coded as 10)
> # otherwise will make paul and cruz look too far right
> newcces <- cces[!is.na(cces$self_emoti) & !is.na(cces$self_consc) & !is.na(cces$self_agree) 
+                 & !is.na(cces$self_extra) & !is.na(cces$self_openn) & !is.na(cces$CC421a) 
+                 & !is.na(cces$newsint) 
+                 & !(cces$CC334C == "Not Asked") 
+                 & !(cces$CC334D == "Not Asked") 
+                 & !(cces$CC334E == "Not Asked") 
+                 & !(cces$CC334F == "Not Asked") 
+                 & !(cces$CC334G == "Not Asked") 
+                 & !(cces$CC334K == "Not Asked") 
+                 & !(cces$CC334L == "Not Asked") 
+                 & !(cces$CC334M == "Not Asked") 
+                 & !(cces$CC334W == "Not Asked"),]
>                 
> # pulling out the stimuli
> aldmckspace <- cbind(as.numeric(newcces$CC334A), # self
+                      as.numeric(newcces$CC334C), # obama
+                      as.numeric(newcces$CC334D), # clinton
+                      as.numeric(newcces$CC334E), # cruz
+                      as.numeric(newcces$CC334F), # paul
+                      as.numeric(newcces$CC334G), # bush
+                      as.numeric(newcces$CC334K), # demparty
+                      as.numeric(newcces$CC334L), # repparty
+                      as.numeric(newcces$CC334M), # teaparty
+                      as.numeric(newcces$CC334W)) # scotus
> 
> # naming the columns/stimuli
> colnames(aldmckspace) <- c("Self", "Obama", "Clinton", "Cruz", "Paul", "Bush", "Democrats", "Republicans", "Tea Party", "Supreme Court")
> 
> # estimating the aldrich-mckelvey model
> aldmck.out <- aldmck(aldmckspace, respondent = 1, polarity = 2, missing = c(8, 9))
> 
> # extracting the estimated "political information" from the A-M results
> newcces$aldmck_polinfo <- aldmck.out$respondents$polinfo # aldrich-mckelvey political information (-1,1 correlation)
> 
> 
> # tobit models with dependent variable as correlation (political information)
> tobit.mod.1 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
+                                       self_agree + self_extra + self_openn, 
+                                       data = newcces, left = -1, right = 1)
>                                       
> tobit.mod.2 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
+                                       self_agree + self_extra + self_openn + 
+                                       I(gender == "Female") + I(2014 - birthyr) + 
+                                       I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                       I(race == "Hispanic") + 
+                                       I(race == "Asian" | race == "Native American" | 
+                                         race == "Mixed" | race == "Other" | race == "Middle Eastern"), 
+                                       data = newcces, left = -1, right = 1)
>                                       
> tobit.mod.3 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
+                                       self_agree + self_extra + self_openn + 
+                                       I(gender == "Female") + I(2014 - birthyr) + 
+                                       I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                       I(race == "Hispanic") + 
+                                       I(race == "Asian" | race == "Native American" | 
+                                         race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                       as.numeric(educ), data = newcces, left = -1, right = 1)
>                                       
> tobit.mod.4 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
+                                       self_agree + self_extra + self_openn + 
+                                       I(gender == "Female") + I(2014 - birthyr) + 
+                                       I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                       I(race == "Hispanic") + 
+                                       I(race == "Asian" | race == "Native American" | 
+                                         race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                       as.numeric(educ) + I(newsint == "Most of the time") + 
+                                       I(newsint == "Don't know"), data = newcces, left = -1, right = 1)
>                                       
> tobit.mod.5 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
+                                       self_agree + self_extra + self_openn + 
+                                       I(gender == "Female") + I(2014 - birthyr) + 
+                                       I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                       I(race == "Hispanic") + 
+                                       I(race == "Asian" | race == "Native American" | 
+                                         race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                       as.numeric(educ) + I(newsint == "Most of the time") + 
+                                       I(newsint == "Don't know") + as.numeric(faminc_recode) + 
+                                       I(faminc_recode == "Prefer not to say"), data = newcces, left = -1, right = 1)
>                                       
> tobit.mod.6 <- tobit(aldmck_polinfo ~ I(1 - self_emoti) + self_consc + 
+                                       self_agree + self_extra + self_openn + 
+                                       I(gender == "Female") + I(2014 - birthyr) + 
+                                       I((2014 - birthyr)^2/100) + I(race == "Black") + 
+                                       I(race == "Hispanic") + 
+                                       I(race == "Asian" | race == "Native American" | 
+                                         race == "Mixed" | race == "Other" | race == "Middle Eastern") + 
+                                       as.numeric(educ) + I(newsint == "Most of the time") + 
+                                       I(newsint == "Don't know") + as.numeric(faminc_recode) + 
+                                       I(faminc_recode == "Prefer not to say") + 
+                                       I(employ == "Full-time") + I(employ == "Part-time") + 
+                                       I(employ == "Unemployed") + I(employ == "Retired"), 
+                                       data = newcces, left = -1, right = 1)
> 
> # coercing the tobit models into a list for texreg
> tobit.mod.list <- list(tobit.mod.1, 
+                        tobit.mod.2,
+                        tobit.mod.3,
+                        tobit.mod.4,
+                        tobit.mod.5,
+                        tobit.mod.6)
> 
> ## table 2
> # texreg for latex (commented out since it's already in the paper)
> # texreg(tobit.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
>                           # custom.coef.names = c("Constant",
>                                                 # "Neuroticism",
>                                                 # "Conscientiousness",
>                                                 # "Agreeableness",
>                                                 # "Extraversion",
>                                                 # "Openness",
>                                                 # "Ln(scale)",
>                                                 # "Female",
>                                                 # "Age",
>                                                 # "Age^2/100",
>                                                 # "Black",
>                                                 # "Hispanic",
>                                                 # "Other Race",
>                                                 # "Education (1 = No HS; 6 = Postgrad)",
>                                                 # "High News Interest",
>                                                 # "Unknown News Interest",
>                                                 # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
>                                                 # "Income Refused",
>                                                 # "Employed Full-Time",
>                                                 # "Employed Part-Time",
>                                                 # "Unemployed",
>                                                 # "Retired"),
>                           # reorder.coef = c(2:6, 8:22,1,7),
>                           # include.aic = FALSE, 
>                           # include.deviance = FALSE)
> 
> # screenreg for review
> # screenreg(tobit.mod.list, digits = 3, stars = c(0.01, 0.05, 0.1),
>                           # custom.coef.names = c("Constant",
>                                                 # "Neuroticism",
>                                                 # "Conscientiousness",
>                                                 # "Agreeableness",
>                                                 # "Extraversion",
>                                                 # "Openness",
>                                                 # "Ln(scale)",
>                                                 # "Female",
>                                                 # "Age",
>                                                 # "Age^2/100",
>                                                 # "Black",
>                                                 # "Hispanic",
>                                                 # "Other Race",
>                                                 # "Education (1 = No HS; 6 = Postgrad)",
>                                                 # "High News Interest",
>                                                 # "Unknown News Interest",
>                                                 # "Income (1 = <10k; 12 = >150k; 13 = Refused)",
>                                                 # "Income Refused",
>                                                 # "Employed Full-Time",
>                                                 # "Employed Part-Time",
>                                                 # "Unemployed",
>                                                 # "Retired"),
>                           # reorder.coef = c(2:6, 8:22,1,7),
>                           # include.aic = FALSE, 
>                           # include.deviance = FALSE,
>                           # include.wald = FALSE)
> 
> 
> # html table
> htmlreg(tobit.mod.list, file = "table2.html", digits = 3, stars = c(0.01, 0.05, 0.1),
+                           custom.coef.names = c("Constant",
+                                                 "Neuroticism",
+                                                 "Conscientiousness",
+                                                 "Agreeableness",
+                                                 "Extraversion",
+                                                 "Openness",
+                                                 "Ln(scale)",
+                                                 "Female",
+                                                 "Age",
+                                                 "Age^2/100",
+                                                 "Black",
+                                                 "Hispanic",
+                                                 "Other Race",
+                                                 "Education (1 = No HS; 6 = Postgrad)",
+                                                 "High News Interest",
+                                                 "Unknown News Interest",
+                                                 "Income (1 = <10k; 12 = >150k; 13 = Refused)",
+                                                 "Income Refused",
+                                                 "Employed Full-Time",
+                                                 "Employed Part-Time",
+                                                 "Unemployed",
+                                                 "Retired"),
+                           reorder.coef = c(2:6, 8:22,1,7),
+                           include.aic = FALSE, 
+                           include.deviance = FALSE,
+                           include.wald = FALSE)
The table was written to the file 'table2.html'.

> 
> 
> # finding the tobit model with the smallest BIC
> tobit.min.bic <- which(c(BIC(tobit.mod.1), BIC(tobit.mod.2), 
+                          BIC(tobit.mod.3), BIC(tobit.mod.4), 
+                          BIC(tobit.mod.5), BIC(tobit.mod.6)) == 
+                        min(c(BIC(tobit.mod.1), BIC(tobit.mod.2), 
+                              BIC(tobit.mod.3), BIC(tobit.mod.4), 
+                              BIC(tobit.mod.5), BIC(tobit.mod.6))))
> 
> # setting up the covariate values for prediction
> tobit.preds <- predict(tobit.mod.list[[tobit.min.bic]], 
+                        newdata = data.frame(self_emoti = seq(mean(newcces$self_emoti) - 2*sd(newcces$self_emoti), 
+                                                              mean(newcces$self_emoti) + 2*sd(newcces$self_emoti), 
+                                                              length.out = 100),
+                                             self_consc = mean(newcces$self_consc),
+                                             self_agree = mean(newcces$self_agree),
+                                             self_extra = mean(newcces$self_extra),
+                                             self_openn = mean(newcces$self_openn),
+                                             newsint = levels(newcces$newsint)[median(as.numeric(newcces$newsint))],
+                                             educ = Mode(newcces$educ),
+                                             employ = Mode(newcces$employ),
+                                             faminc_recode = Mode(newcces$faminc_recode),
+                                             gender = Mode(newcces$gender),
+                                             birthyr = median(newcces$birthyr),
+                                             race = Mode(newcces$race)), se = TRUE, type = "linear")
> 
> # estimating the predicted correlations with 90% CIs
> # reversed scale because moving from emotional stability to neuroticism
> tobit.predvals <- data.frame(neuroticism = seq(2, -2, length.out = 100),
+                                 low = tobit.preds$fit - qnorm(0.95)*tobit.preds$se,
+                                 mid = tobit.preds$fit,
+                                 high = tobit.preds$fit + qnorm(0.95)*tobit.preds$se)
> 
> 
> # combining the tobit plot with the binomial plot from before
> binom.tobit.predvals <- rbind(data.frame(frame = total.na.predvals, plot = 1),
+                               data.frame(frame = tobit.predvals, plot = 2))
> 
> # naming the columns of the data frame
> colnames(binom.tobit.predvals) <- c("Neuroticism",
+                                     "low",
+                                     "mid",
+                                     "high",
+                                     "DV")
> 
> # setting up the appropriate result labels
> binom.tobit.predvals$DV <- factor(binom.tobit.predvals$DV, 
+                                   labels = c("Proportion of NA/DK Responses",
+                                              "Correlation of Perceived Ideological Space with True Space"))
> 
> # figure 4
> combined.plot <- ggplot(binom.tobit.predvals, aes(x = Neuroticism, y = mid, linetype = DV)) + 
+                     geom_line(size = 0.75) + 
+                     geom_ribbon(aes(ymin = low, ymax = high, fill = DV, linetype = NA), alpha = 0.6) + 
+                     theme_bw() + 
+                     xlab("Neuroticism\n(Number of Standard Deviations From Mean)") + 
+                     ylab("Predicted Value of Dependent Variable\n") + 
+                     scale_colour_discrete(name = c("Dependent Variable")) + 
+                     scale_fill_grey(name = c("Dependent Variable"), start = 0.4, end = 0.8) +
+                     scale_linetype_discrete(name = c("Dependent Variable")) + 
+                     ylim(0,1) + 
+                     theme(legend.position="bottom", legend.direction = "vertical")
> 
> 
> cairo_ps(file="combined_plot.eps", height = 5, width = 7.5)
> combined.plot
> dev.off()
null device 
          1 
> 
> 
> 
> 
> ### load hierarchical model data
> load("DKWYGCore1.Rda") 
> load("DKWYGCore2.Rda") 
> load("DKWYGCore3.Rda") 
> load("DKWYGCore4.Rda") 
> load("DKWYGCore5.Rda") 
> load("DKWYGCore6.Rda")
> load("DKWYGCore7.Rda") 
> load("DKWYGCore8.Rda") 
> 
> thinning <- 4 # 12,500 per chain; 100,000 overall
> a <- seq(1, dim(dkwyg.1[[1]])[1])
> b <- a[seq(1, length(a), thinning)]
> 
> 
> # coerce into a single object
> DKWYG <- mcmc.list(
+                    as.mcmc(dkwyg.1[[1]][b,]),
+                    as.mcmc(dkwyg.2[[1]][b,]),
+                    as.mcmc(dkwyg.3[[1]][b,]),
+                    as.mcmc(dkwyg.4[[1]][b,]),
+                    as.mcmc(dkwyg.5[[1]][b,]),
+                    as.mcmc(dkwyg.6[[1]][b,]),
+                    as.mcmc(dkwyg.7[[1]][b,]),
+                    as.mcmc(dkwyg.8[[1]][b,])
+                    )
> 
> 
> 
> # coerce the chains into a single chain for further analysis
> BZ_data <- as.data.frame(runjags::combine.mcmc(DKWYG))
> 
> candNames <- c("Obama","Cruz","Clinton","Paul","Bush","Democrats",
+                "Republicans","Tea Party","Supreme Court")
> colnames(BZ_data)[86:94] <- candNames
> melted_BZ <- melt(BZ_data)
No id variables; using all as measure variables
> 
> # calculating posterior probability of parameter having same sign as posterior mean
> postprob <- function(x){
+    foo <- sum(sign(x) == sign(mean(x)))/length(x)
+    return(foo)
+    }
> 
> round(data.frame(apply(BZ_data, 2, postprob)), digits = 3)
              apply.BZ_data..2..postprob.
beta_a[1]                           0.781
beta_a[2]                           1.000
beta_a[3]                           0.849
beta_a[4]                           0.986
beta_a[5]                           0.923
beta_a[6]                           0.970
beta_a[7]                           0.594
beta_a[8]                           0.558
beta_a[9]                           0.552
beta_a[10]                          1.000
beta_a[11]                          0.875
beta_a[12]                          0.936
beta_a[13]                          0.995
beta_a[14]                          0.723
beta_a[15]                          0.799
beta_a[16]                          0.913
beta_a[17]                          0.552
beta_a[18]                          0.628
beta_a[19]                          0.964
beta_a[20]                          0.545
beta_a[21]                          0.559
beta_b[1]                           0.976
beta_b[2]                           0.580
beta_b[3]                           0.998
beta_b[4]                           0.802
beta_b[5]                           0.873
beta_b[6]                           0.984
beta_b[7]                           0.708
beta_b[8]                           0.933
beta_b[9]                           0.968
beta_b[10]                          0.989
beta_b[11]                          0.651
beta_b[12]                          0.513
beta_b[13]                          0.999
beta_b[14]                          1.000
beta_b[15]                          0.788
beta_b[16]                          0.510
beta_b[17]                          0.889
beta_b[18]                          0.729
beta_b[19]                          0.646
beta_b[20]                          0.995
beta_b[21]                          0.774
beta_d[1]                           0.905
beta_d[2]                           0.989
beta_d[3]                           0.962
beta_d[4]                           0.996
beta_d[5]                           0.816
beta_d[6]                           0.712
beta_d[7]                           0.963
beta_d[8]                           0.791
beta_d[9]                           0.859
beta_d[10]                          1.000
beta_d[11]                          0.999
beta_d[12]                          0.854
beta_d[13]                          0.974
beta_d[14]                          1.000
beta_d[15]                          0.772
beta_d[16]                          0.736
beta_d[17]                          0.967
beta_d[18]                          0.831
beta_d[19]                          0.592
beta_d[20]                          0.737
beta_d[21]                          0.692
beta_e[1]                           0.790
beta_e[2]                           0.609
beta_e[3]                           0.847
beta_e[4]                           0.500
beta_e[5]                           0.965
beta_e[6]                           0.855
beta_e[7]                           1.000
beta_e[8]                           0.893
beta_e[9]                           0.675
beta_e[10]                          0.833
beta_e[11]                          0.946
beta_e[12]                          0.971
beta_e[13]                          1.000
beta_e[14]                          1.000
beta_e[15]                          1.000
beta_e[16]                          0.989
beta_e[17]                          0.771
beta_e[18]                          0.803
beta_e[19]                          0.654
beta_e[20]                          0.896
beta_e[21]                          0.855
beta_s                              1.000
Obama                               1.000
Cruz                                1.000
Clinton                             1.000
Paul                                1.000
Bush                                1.000
Democrats                           1.000
Republicans                         1.000
Tea Party                           1.000
Supreme Court                       0.548
tau[1]                              0.664
tau[2]                              1.000
tau[3]                              1.000
tau[4]                              1.000
tau[5]                              1.000
> 
>    
> # functions to create HPDs
> middle95HPD <- function(x){
+     out <- HPDinterval(as.mcmc(x), prob = 0.95)[,]
+     names(out) <- c("ymin","ymax")
+     return(out)
+ }
> 
> # alternative with middle 80%
> middle80HPD <- function(x){
+     out <- HPDinterval(as.mcmc(x), prob = 0.8)[,]
+     names(out) <- c("ymin","ymax")
+     return(out)
+ }
> 
> # alternative with middle 75%
> middle75HPD <- function(x){
+     out <- HPDinterval(as.mcmc(x), prob = 0.75)[,]
+     names(out) <- c("ymin","ymax")
+     return(out)
+ }
> 
> 
> # alternative with middle 50%
> middle50HPD <- function(x){
+     out <- HPDinterval(as.mcmc(x), prob = 0.5)[,]
+     names(out) <- c("ymin","ymax")
+     return(out)
+ }
> 
> 
> 
> ### distributional plots
> # load bootstrapped AM replications
> load("AMBoot.Rda") 
> 
> # rescale the replications so Obama = -1 and Cruz = 1
> aldmck.boot.rescale <- 2*(aldmck.boot - aldmck.boot[,1])/(aldmck.boot[,3] - aldmck.boot[,1]) - 1
> 
> # use first 100,000 points since the thinned chains combined have 100,000 draws
> aldmck.boot.rescale <- aldmck.boot.rescale[1:100000,]
> 
> 
> # big plot for comparing distributions of estimates
> melt.aldmck <- melt(aldmck.boot.rescale[,])
> melt.aldmck <- melt.aldmck[,-1]
> colnames(melt.aldmck)[1] <- "variable"
> melt.aldmck$variable <- reorder(melt.aldmck$variable, melt.aldmck$value)
> melted_BZ$variable <- reorder(melted_BZ$variable, melted_BZ$value)
> 
> # make sure order is correct!
> aldmck.means <- aggregate(melt.aldmck[,2], list(melt.aldmck[,1]), mean)
> aldmck.means <- aldmck.means[aldmck.means$Group.1 %in% candNames[-c(1,2)],]
> aldmck.means <- aldmck.means[order(as.character(aldmck.means$Group.1)),]
> BZ.means <- aggregate(melted_BZ[which(melted_BZ$variable %in% candNames),2], 
+                       list(melted_BZ[which(melted_BZ$variable %in% candNames),1]), mean) 
> BZ.means <- BZ.means[BZ.means$Group.1 %in% candNames[-c(1,2)],]
> BZ.means <- BZ.means[order(as.character(BZ.means$Group.1)),]
> r.squared <- summary(lm(aldmck.means[,2] ~ BZ.means[,2]))$r.squared
> 
> # combining the AM and BZ estimates into a data frame and labeling the methods accordingly
> est_comparison_data <- rbind(data.frame(melt.aldmck, Method = "Aldrich-McKelvey"),
+                              data.frame(melted_BZ[melted_BZ$variable %in% candNames,], Method = "Hierarchical Model"))
> 
> # figure 5
> big_comparison_plot <- ggplot(est_comparison_data, aes(x=variable, y=value, fill = Method, shape = Method)) + 
+                           stat_summary(data = est_comparison_data[est_comparison_data$variable %in% candNames[-c(1:2)],],
+                                        fun.y=median,geom='point',position=position_dodge(width=0.725)) +
+                           geom_violin(data = est_comparison_data[est_comparison_data$variable %in% candNames[-c(1:2)],], 
+                                       aes(x= variable,y=value,fill=Method), alpha = 0.6, width = 0.75) + 
+                           stat_summary(data = est_comparison_data[est_comparison_data$variable %in% candNames[-c(1:2)],],
+                                        fun.y=median,geom='point',position=position_dodge(width=0.725)) +
+                           stat_summary(data = est_comparison_data[est_comparison_data$variable %in% candNames[c(1:2)],],
+                                        fun.y=median,geom='point') +
+                           coord_flip() + 
+                           theme_bw() +
+                           theme(axis.title.y = element_blank()) +
+                           scale_x_discrete(limits = rev(levels(factor(melt.aldmck[which(melt.aldmck$variable %in% candNames),]$variable)))) + 
+                           scale_fill_manual(values = c("grey40", "white"), name = "Method") + 
+                           ylab("\nEstimated Ideology") + 
+                           annotate("text", x = 8, y = 1, label = paste("R^2 %~~%", sprintf("%.3f",round(r.squared, digits = 4))), parse = TRUE)
> 
> 
> cairo_ps(file="bigcomparisonplot.eps", height = 5, width = 7.5)
> big_comparison_plot
> dev.off()
null device 
          1 
> 
> 
> 
> 
> # examining the results of the hierarchical model
> personality <- c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", 
+               "Neuroticism")
> personalityNews <- c("Openness", "Conscientiousness", "Extraversion", "Agreeableness", 
+               "Neuroticism","High News Interest")
> varNames <- c("Constant", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", 
+               "Neuroticism", "Female",
+               "Age", "Age^2/100", "Black", "Hispanic", "Other Race",
+               "Education", "High News Interest", "Unknown News Interest",
+               "Income", "Income Refused",
+               "Full Time Employment", "Part Time Employment", "Unemployed", "Retired")
> 
> colnames(BZ_data)[1:(length(varNames))] <- paste("beta_a[",varNames,"]",sep="")
> colnames(BZ_data)[(length(varNames)+1):(2*length(varNames))] <- paste("beta_b[",varNames,"]",sep="")
> colnames(BZ_data)[(2*length(varNames)+1):(3*length(varNames))] <- paste("beta_d[",varNames,"]",sep="")
> colnames(BZ_data)[(3*length(varNames)+1):(4*length(varNames))] <- paste("beta_e[",varNames,"]",sep="")
> 
> 
> # coerce into a separate object for testing
> DKWYG.nofixed <- mcmc.list(
+                            as.mcmc(dkwyg.1[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.2[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.3[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.4[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.5[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.6[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.7[[1]][b,-c(86,87)]),
+                            as.mcmc(dkwyg.8[[1]][b,-c(86,87)])
+                            )
> 
> gelman.diag(DKWYG.nofixed)
Potential scale reduction factors:

           Point est. Upper C.I.
beta_a[1]        1.00       1.00
beta_a[2]        1.00       1.00
beta_a[3]        1.00       1.00
beta_a[4]        1.00       1.00
beta_a[5]        1.00       1.00
beta_a[6]        1.00       1.00
beta_a[7]        1.00       1.00
beta_a[8]        1.00       1.00
beta_a[9]        1.00       1.00
beta_a[10]       1.00       1.00
beta_a[11]       1.00       1.00
beta_a[12]       1.00       1.00
beta_a[13]       1.00       1.00
beta_a[14]       1.00       1.00
beta_a[15]       1.00       1.01
beta_a[16]       1.00       1.00
beta_a[17]       1.00       1.00
beta_a[18]       1.00       1.00
beta_a[19]       1.00       1.00
beta_a[20]       1.00       1.00
beta_a[21]       1.00       1.00
beta_b[1]        1.00       1.00
beta_b[2]        1.00       1.00
beta_b[3]        1.00       1.00
beta_b[4]        1.00       1.00
beta_b[5]        1.00       1.00
beta_b[6]        1.00       1.00
beta_b[7]        1.00       1.00
beta_b[8]        1.00       1.00
beta_b[9]        1.00       1.00
beta_b[10]       1.00       1.00
beta_b[11]       1.00       1.00
beta_b[12]       1.00       1.00
beta_b[13]       1.00       1.00
beta_b[14]       1.00       1.00
beta_b[15]       1.00       1.01
beta_b[16]       1.00       1.00
beta_b[17]       1.00       1.00
beta_b[18]       1.00       1.00
beta_b[19]       1.00       1.00
beta_b[20]       1.00       1.00
beta_b[21]       1.00       1.00
beta_d[1]        1.00       1.00
beta_d[2]        1.02       1.04
beta_d[3]        1.01       1.03
beta_d[4]        1.00       1.01
beta_d[5]        1.01       1.01
beta_d[6]        1.01       1.03
beta_d[7]        1.01       1.02
beta_d[8]        1.01       1.02
beta_d[9]        1.01       1.01
beta_d[10]       1.01       1.01
beta_d[11]       1.00       1.01
beta_d[12]       1.03       1.06
beta_d[13]       1.01       1.03
beta_d[14]       1.01       1.03
beta_d[15]       1.01       1.02
beta_d[16]       1.02       1.03
beta_d[17]       1.02       1.05
beta_d[18]       1.01       1.03
beta_d[19]       1.01       1.02
beta_d[20]       1.01       1.02
beta_d[21]       1.03       1.06
beta_e[1]        1.00       1.00
beta_e[2]        1.01       1.01
beta_e[3]        1.01       1.01
beta_e[4]        1.00       1.01
beta_e[5]        1.01       1.02
beta_e[6]        1.00       1.01
beta_e[7]        1.01       1.03
beta_e[8]        1.01       1.01
beta_e[9]        1.01       1.01
beta_e[10]       1.02       1.04
beta_e[11]       1.01       1.02
beta_e[12]       1.01       1.01
beta_e[13]       1.02       1.04
beta_e[14]       1.02       1.04
beta_e[15]       1.01       1.01
beta_e[16]       1.02       1.03
beta_e[17]       1.02       1.04
beta_e[18]       1.00       1.01
beta_e[19]       1.00       1.01
beta_e[20]       1.00       1.01
beta_e[21]       1.02       1.05
beta_s           1.04       1.09
gamma[3]         1.00       1.00
gamma[4]         1.00       1.01
gamma[5]         1.01       1.02
gamma[6]         1.00       1.00
gamma[7]         1.01       1.02
gamma[8]         1.00       1.01
gamma[9]         1.01       1.01
tau[1]           1.01       1.02
tau[2]           1.01       1.02
tau[3]           1.01       1.03
tau[4]           1.02       1.04
tau[5]           1.01       1.03

Multivariate psrf

1.1
> 
> # print table a-1
> table_a1_out <- gelman.diag(DKWYG.nofixed)
> table_a1_psrf <- data.frame(table_a1_out$psrf[,2])
> table_a1 <- data.frame(c(table_a1_psrf[[1]], table_a1_out$mpsrf))
> rownames(table_a1) <- c(rownames(table_a1_psrf), "Multivariate PSRF")
> rownames(table_a1)[1:(length(varNames))] <- paste("beta_a[",varNames,"]",sep="")
> rownames(table_a1)[(length(varNames)+1):(2*length(varNames))] <- paste("beta_b[",varNames,"]",sep="")
> rownames(table_a1)[(2*length(varNames)+1):(3*length(varNames))] <- paste("beta_d[",varNames,"]",sep="")
> rownames(table_a1)[(3*length(varNames)+1):(4*length(varNames))] <- paste("beta_e[",varNames,"]",sep="")
> rownames(table_a1) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a1)))))))
> rownames(table_a1)[rownames(table_a1) == "c[5]"] <- "c[6]"
> rownames(table_a1)[rownames(table_a1) == "c[4]"] <- "c[5]"
> rownames(table_a1)[rownames(table_a1) == "c[3]"] <- "c[4]"
> rownames(table_a1)[rownames(table_a1) == "c[2]"] <- "c[3]"
> rownames(table_a1)[rownames(table_a1) == "c[1]"] <- "c[2]"
> rownames(table_a1)[rownames(table_a1) == "gamma[3]"] <- "Clinton"
> rownames(table_a1)[rownames(table_a1) == "gamma[4]"] <- "Paul"
> rownames(table_a1)[rownames(table_a1) == "gamma[5]"] <- "Bush"
> rownames(table_a1)[rownames(table_a1) == "gamma[6]"] <- "Democrats"
> rownames(table_a1)[rownames(table_a1) == "gamma[7]"] <- "Republicans"
> rownames(table_a1)[rownames(table_a1) == "gamma[8]"] <- "Tea Party"
> rownames(table_a1)[rownames(table_a1) == "gamma[9]"] <- "Supreme Court"
> colnames(table_a1) <- "PSRF Upper Bound"
> 
> print(xtable(table_a1), type = "html", file = "table-a1.html")
> 
> 
> 
> # hpd intervals
> BZ_sum <- summary(as.mcmc(BZ_data))
> hpd99 <- HPDinterval(as.mcmc(BZ_data), prob = 0.99)
> hpd95 <- HPDinterval(as.mcmc(BZ_data), prob = 0.95)
> hpd90 <- HPDinterval(as.mcmc(BZ_data), prob = 0.9)
> hpd80 <- HPDinterval(as.mcmc(BZ_data), prob = 0.8)
> hpd70 <- HPDinterval(as.mcmc(BZ_data), prob = 0.7)
> hpd60 <- HPDinterval(as.mcmc(BZ_data), prob = 0.6)
> hpd50 <- HPDinterval(as.mcmc(BZ_data), prob = 0.5)
> 
> hpd99.sig <- as.numeric(apply(hpd99,1, function(x){0 > max(x) | 0 < min(x)}))
> hpd95.sig <- as.numeric(apply(hpd95,1, function(x){0 > max(x) | 0 < min(x)}))
> hpd90.sig <- as.numeric(apply(hpd90,1, function(x){0 > max(x) | 0 < min(x)}))
> hpd80.sig <- as.numeric(apply(hpd80,1, function(x){0 > max(x) | 0 < min(x)}))
> hpd70.sig <- as.numeric(apply(hpd70,1, function(x){0 > max(x) | 0 < min(x)}))
> hpd60.sig <- as.numeric(apply(hpd60,1, function(x){0 > max(x) | 0 < min(x)}))
> hpd50.sig <- as.numeric(apply(hpd50,1, function(x){0 > max(x) | 0 < min(x)}))
> 
> smallestHPD <- apply(cbind(hpd99.sig, hpd95.sig, hpd90.sig, hpd80.sig, hpd70.sig, hpd60.sig, hpd50.sig),1,which.max)
> neversig <- as.numeric(apply(cbind(hpd99.sig, hpd95.sig, hpd90.sig, hpd80.sig, hpd70.sig, hpd60.sig, hpd50.sig),1,sum) != 0)
> smallestHPD[smallestHPD*neversig == 0] <- NA
> smallestHPD <- as.numeric(gsub(".sig","",gsub("hpd","",colnames(cbind(hpd99.sig, hpd95.sig, hpd90.sig, hpd80.sig, hpd70.sig, hpd60.sig, hpd50.sig))[smallestHPD])))
> smallestHPD.result <- data.frame(mean = BZ_sum$statistics[,'Mean'], se = BZ_sum$statistics[,'SD'], hpd80low = hpd80[,1], hpd80high = hpd80[,2], smallestHPD = smallestHPD)
> rownames(smallestHPD.result) <- gsub("beta_s","sal_slo",gsub("beta_e","sal_int",gsub("beta_d","decisive",gsub("beta_b","op_slo",gsub("beta_a","op_int",rownames(smallestHPD.result))))))
> 
> # hpd intervals mentioned in text
> round(HPDinterval(as.mcmc(BZ_data), prob = 0.69), digits = 3)
                               lower  upper
beta_a[Constant]              -0.170  1.298
beta_a[Openness]               1.394  2.140
beta_a[Conscientiousness]     -0.763 -0.017
beta_a[Extraversion]          -0.914 -0.336
beta_a[Agreeableness]          0.165  0.950
beta_a[Neuroticism]            0.278  0.945
beta_a[Female]                -0.175  0.106
beta_a[Age]                   -0.021  0.028
beta_a[Age^2/100]             -0.029  0.023
beta_a[Black]                  0.665  1.112
beta_a[Hispanic]               0.041  0.586
beta_a[Other Race]            -0.619 -0.122
beta_a[Education]              0.074  0.172
beta_a[High News Interest]    -0.233  0.062
beta_a[Unknown News Interest] -0.999  0.108
beta_a[Income]                -0.059 -0.009
beta_a[Income Refused]        -0.228  0.297
beta_a[Full Time Employment]  -0.128  0.249
beta_a[Part Time Employment]  -0.667 -0.188
beta_a[Unemployed]            -0.359  0.275
beta_a[Retired]               -0.284  0.205
beta_b[Constant]               0.729  2.260
beta_b[Openness]              -0.321  0.461
beta_b[Conscientiousness]      0.723  1.509
beta_b[Extraversion]          -0.552  0.058
beta_b[Agreeableness]          0.054  0.880
beta_b[Neuroticism]           -1.103 -0.401
beta_b[Female]                -0.228  0.070
beta_b[Age]                   -0.065 -0.013
beta_b[Age^2/100]              0.024  0.078
beta_b[Black]                 -0.783 -0.310
beta_b[Hispanic]              -0.396  0.187
beta_b[Other Race]            -0.247  0.271
beta_b[Education]              0.115  0.218
beta_b[High News Interest]     0.544  0.860
beta_b[Unknown News Interest] -1.020  0.146
beta_b[Income]                -0.026  0.027
beta_b[Income Refused]         0.063  0.616
beta_b[Full Time Employment]  -0.077  0.320
beta_b[Part Time Employment]  -0.152  0.353
beta_b[Unemployed]             0.539  1.209
beta_b[Retired]               -0.448  0.067
beta_d[Constant]               0.339  2.891
beta_d[Openness]               0.943  2.457
beta_d[Conscientiousness]     -2.563 -0.708
beta_d[Extraversion]          -2.071 -0.928
beta_d[Agreeableness]         -0.108  1.503
beta_d[Neuroticism]           -0.956  0.294
beta_d[Female]                -0.881 -0.232
beta_d[Age]                   -0.089  0.012
beta_d[Age^2/100]              0.001  0.106
beta_d[Black]                  1.465  3.460
beta_d[Hispanic]               1.416  3.655
beta_d[Other Race]            -0.270  1.519
beta_d[Education]              0.088  0.302
beta_d[High News Interest]     1.028  1.701
beta_d[Unknown News Interest] -0.549  2.164
beta_d[Income]                -0.024  0.093
beta_d[Income Refused]        -1.660 -0.482
beta_d[Full Time Employment]  -0.020  0.673
beta_d[Part Time Employment]  -0.362  0.531
beta_d[Unemployed]            -0.417  1.098
beta_d[Retired]               -0.267  0.736
beta_e[Constant]              -2.309  0.305
beta_e[Openness]              -0.615  1.079
beta_e[Conscientiousness]     -0.008  1.610
beta_e[Extraversion]          -0.695  0.648
beta_e[Agreeableness]         -2.501 -0.736
beta_e[Neuroticism]           -1.445 -0.027
beta_e[Female]                -1.485 -0.806
beta_e[Age]                    0.012  0.114
beta_e[Age^2/100]             -0.082  0.030
beta_e[Black]                 -0.891  0.042
beta_e[Hispanic]              -1.311 -0.278
beta_e[Other Race]            -1.537 -0.470
beta_e[Education]              0.318  0.574
beta_e[High News Interest]     1.773  2.532
beta_e[Unknown News Interest] -3.330 -1.921
beta_e[Income]                 0.070  0.194
beta_e[Income Refused]        -1.267  0.188
beta_e[Full Time Employment]  -0.765  0.076
beta_e[Part Time Employment]  -0.321  0.710
beta_e[Unemployed]            -1.292 -0.138
beta_e[Retired]               -1.225 -0.013
beta_s                         0.400  0.726
Obama                         -1.000 -1.000
Cruz                           1.000  1.000
Clinton                       -0.935 -0.872
Paul                           0.803  0.871
Bush                           0.563  0.625
Democrats                     -0.991 -0.926
Republicans                    0.642  0.706
Tea Party                      1.251  1.340
Supreme Court                 -0.030  0.024
tau[1]                        -0.084  0.033
tau[2]                         1.197  1.343
tau[3]                         2.503  2.673
tau[4]                         3.511  3.697
tau[5]                         5.256  5.472
attr(,"Probability")
[1] 0.69
> round(HPDinterval(as.mcmc(BZ_data), prob = 0.41), digits = 3)
                               lower  upper
beta_a[Constant]               0.205  0.986
beta_a[Openness]               1.578  1.974
beta_a[Conscientiousness]     -0.584 -0.186
beta_a[Extraversion]          -0.791 -0.485
beta_a[Agreeableness]          0.345  0.762
beta_a[Neuroticism]            0.437  0.791
beta_a[Female]                -0.103  0.046
beta_a[Age]                   -0.011  0.015
beta_a[Age^2/100]             -0.016  0.011
beta_a[Black]                  0.777  1.015
beta_a[Hispanic]               0.156  0.445
beta_a[Other Race]            -0.507 -0.244
beta_a[Education]              0.097  0.149
beta_a[High News Interest]    -0.166 -0.009
beta_a[Unknown News Interest] -0.775 -0.188
beta_a[Income]                -0.046 -0.020
beta_a[Income Refused]        -0.104  0.173
beta_a[Full Time Employment]  -0.040  0.158
beta_a[Part Time Employment]  -0.566 -0.312
beta_a[Unemployed]            -0.196  0.139
beta_a[Retired]               -0.166  0.094
beta_b[Constant]               1.093  1.904
beta_b[Openness]              -0.147  0.269
beta_b[Conscientiousness]      0.909  1.323
beta_b[Extraversion]          -0.423 -0.102
beta_b[Agreeableness]          0.247  0.684
beta_b[Neuroticism]           -0.938 -0.565
beta_b[Female]                -0.162 -0.004
beta_b[Age]                   -0.052 -0.025
beta_b[Age^2/100]              0.036  0.065
beta_b[Black]                 -0.674 -0.424
beta_b[Hispanic]              -0.263  0.046
beta_b[Other Race]            -0.126  0.149
beta_b[Education]              0.137  0.192
beta_b[High News Interest]     0.618  0.785
beta_b[Unknown News Interest] -0.755 -0.141
beta_b[Income]                -0.013  0.015
beta_b[Income Refused]         0.180  0.475
beta_b[Full Time Employment]   0.009  0.221
beta_b[Part Time Employment]  -0.038  0.228
beta_b[Unemployed]             0.695  1.052
beta_b[Retired]               -0.343 -0.069
beta_d[Constant]               0.821  2.172
beta_d[Openness]               1.302  2.097
beta_d[Conscientiousness]     -2.081 -1.093
beta_d[Extraversion]          -1.814 -1.209
beta_d[Agreeableness]          0.242  1.094
beta_d[Neuroticism]           -0.661 -0.003
beta_d[Female]                -0.698 -0.361
beta_d[Age]                   -0.062 -0.008
beta_d[Age^2/100]              0.026  0.082
beta_d[Black]                  1.839  2.878
beta_d[Hispanic]               1.875  3.059
beta_d[Other Race]             0.002  0.883
beta_d[Education]              0.135  0.246
beta_d[High News Interest]     1.173  1.526
beta_d[Unknown News Interest] -0.017  1.411
beta_d[Income]                 0.004  0.065
beta_d[Income Refused]        -1.383 -0.764
beta_d[Full Time Employment]   0.128  0.492
beta_d[Part Time Employment]  -0.168  0.303
beta_d[Unemployed]            -0.058  0.728
beta_d[Retired]               -0.048  0.476
beta_e[Constant]              -1.739 -0.346
beta_e[Openness]              -0.204  0.693
beta_e[Conscientiousness]      0.386  1.245
beta_e[Extraversion]          -0.374  0.332
beta_e[Agreeableness]         -2.059 -1.122
beta_e[Neuroticism]           -1.123 -0.376
beta_e[Female]                -1.319 -0.962
beta_e[Age]                    0.037  0.091
beta_e[Age^2/100]             -0.056  0.004
beta_e[Black]                 -0.653 -0.162
beta_e[Hispanic]              -1.071 -0.524
beta_e[Other Race]            -1.270 -0.703
beta_e[Education]              0.374  0.509
beta_e[High News Interest]     1.953  2.357
beta_e[Unknown News Interest] -2.998 -2.251
beta_e[Income]                 0.100  0.165
beta_e[Income Refused]        -0.952 -0.183
beta_e[Full Time Employment]  -0.564 -0.121
beta_e[Part Time Employment]  -0.091  0.456
beta_e[Unemployed]            -1.033 -0.422
beta_e[Retired]               -0.924 -0.283
beta_s                         0.475  0.648
Obama                         -1.000 -1.000
Cruz                           1.000  1.000
Clinton                       -0.920 -0.887
Paul                           0.817  0.853
Bush                           0.580  0.612
Democrats                     -0.976 -0.941
Republicans                    0.656  0.690
Tea Party                      1.270  1.317
Supreme Court                 -0.019  0.009
tau[1]                        -0.054  0.008
tau[2]                         1.229  1.307
tau[3]                         2.544  2.635
tau[4]                         3.574  3.674
tau[5]                         5.302  5.417
attr(,"Probability")
[1] 0.41
> 
> 
> # tables a-2 through a-8
> table_a2 <- data.frame(hpd99)
> rownames(table_a2) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a2)))))))
> colnames(table_a2) <- c("Lower Bound", "Upper Bound")
> rownames(table_a2)[rownames(table_a2) == "c[5]"] <- "c[6]"
> rownames(table_a2)[rownames(table_a2) == "c[4]"] <- "c[5]"
> rownames(table_a2)[rownames(table_a2) == "c[3]"] <- "c[4]"
> rownames(table_a2)[rownames(table_a2) == "c[2]"] <- "c[3]"
> rownames(table_a2)[rownames(table_a2) == "c[1]"] <- "c[2]"
> print(xtable(table_a2, digits = 3), type = "html", file = "table-a2.html")
> 
> table_a3 <- data.frame(hpd95)
> rownames(table_a3) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a3)))))))
> colnames(table_a3) <- c("Lower Bound", "Upper Bound")
> rownames(table_a3)[rownames(table_a3) == "c[5]"] <- "c[6]"
> rownames(table_a3)[rownames(table_a3) == "c[4]"] <- "c[5]"
> rownames(table_a3)[rownames(table_a3) == "c[3]"] <- "c[4]"
> rownames(table_a3)[rownames(table_a3) == "c[2]"] <- "c[3]"
> rownames(table_a3)[rownames(table_a3) == "c[1]"] <- "c[2]"
> print(xtable(table_a3, digits = 3), type = "html", file = "table-a3.html")
> 
> table_a4 <- data.frame(hpd90)
> rownames(table_a4) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a4)))))))
> colnames(table_a4) <- c("Lower Bound", "Upper Bound")
> rownames(table_a4)[rownames(table_a4) == "c[5]"] <- "c[6]"
> rownames(table_a4)[rownames(table_a4) == "c[4]"] <- "c[5]"
> rownames(table_a4)[rownames(table_a4) == "c[3]"] <- "c[4]"
> rownames(table_a4)[rownames(table_a4) == "c[2]"] <- "c[3]"
> rownames(table_a4)[rownames(table_a4) == "c[1]"] <- "c[2]"
> print(xtable(table_a4, digits = 3), type = "html", file = "table-a4.html")
> 
> table_a5 <- data.frame(hpd80)
> rownames(table_a5) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a5)))))))
> colnames(table_a5) <- c("Lower Bound", "Upper Bound")
> rownames(table_a5)[rownames(table_a5) == "c[5]"] <- "c[6]"
> rownames(table_a5)[rownames(table_a5) == "c[4]"] <- "c[5]"
> rownames(table_a5)[rownames(table_a5) == "c[3]"] <- "c[4]"
> rownames(table_a5)[rownames(table_a5) == "c[2]"] <- "c[3]"
> rownames(table_a5)[rownames(table_a5) == "c[1]"] <- "c[2]"
> print(xtable(table_a5, digits = 3), type = "html", file = "table-a5.html")
> 
> table_a6 <- data.frame(hpd70)
> rownames(table_a6) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a6)))))))
> colnames(table_a6) <- c("Lower Bound", "Upper Bound")
> rownames(table_a6)[rownames(table_a6) == "c[5]"] <- "c[6]"
> rownames(table_a6)[rownames(table_a6) == "c[4]"] <- "c[5]"
> rownames(table_a6)[rownames(table_a6) == "c[3]"] <- "c[4]"
> rownames(table_a6)[rownames(table_a6) == "c[2]"] <- "c[3]"
> rownames(table_a6)[rownames(table_a6) == "c[1]"] <- "c[2]"
> print(xtable(table_a6, digits = 3), type = "html", file = "table-a6.html")
> 
> table_a7 <- data.frame(hpd60)
> rownames(table_a7) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a7)))))))
> colnames(table_a7) <- c("Lower Bound", "Upper Bound")
> rownames(table_a7)[rownames(table_a7) == "c[5]"] <- "c[6]"
> rownames(table_a7)[rownames(table_a7) == "c[4]"] <- "c[5]"
> rownames(table_a7)[rownames(table_a7) == "c[3]"] <- "c[4]"
> rownames(table_a7)[rownames(table_a7) == "c[2]"] <- "c[3]"
> rownames(table_a7)[rownames(table_a7) == "c[1]"] <- "c[2]"
> print(xtable(table_a7, digits = 3), type = "html", file = "table-a7.html")
> 
> table_a8 <- data.frame(hpd50)
> rownames(table_a8) <- gsub("tau", "c", gsub("beta_s","Saliency Slope", gsub("beta_e","Saliency Intercept", gsub("beta_d","Decisiveness",gsub("beta_b","Opinion Slope",gsub("beta_a","Opinion Intercept",rownames(table_a8)))))))
> colnames(table_a8) <- c("Lower Bound", "Upper Bound")
> rownames(table_a8)[rownames(table_a8) == "c[5]"] <- "c[6]"
> rownames(table_a8)[rownames(table_a8) == "c[4]"] <- "c[5]"
> rownames(table_a8)[rownames(table_a8) == "c[3]"] <- "c[4]"
> rownames(table_a8)[rownames(table_a8) == "c[2]"] <- "c[3]"
> rownames(table_a8)[rownames(table_a8) == "c[1]"] <- "c[2]"
> print(xtable(table_a8, digits = 3), type = "html", file = "table-a8.html")
> 
> 
> 
> 
> ### table 3 information 
> table3 <- as.data.frame(smallestHPD.result[!(rownames(smallestHPD.result) %in% candNames),-c(3,4)])
> rownames(table3) <- gsub("tau", "c", gsub("sal_slo","Saliency Slope",gsub("sal_int","Saliency Intercept",gsub("decisive","Decisiveness",gsub("op_slo","Opinion Slope",gsub("op_int","Opinion Intercept",rownames(table3)))))))
> rownames(table3)[rownames(table3) == "c[5]"] <- "c[6]"
> rownames(table3)[rownames(table3) == "c[4]"] <- "c[5]"
> rownames(table3)[rownames(table3) == "c[3]"] <- "c[4]"
> rownames(table3)[rownames(table3) == "c[2]"] <- "c[3]"
> rownames(table3)[rownames(table3) == "c[1]"] <- "c[2]"
> table3 <- rbind(table3, NA) 
> rownames(table3)[dim(table3)[1]] <- "Number of Observations"
> table3[dim(table3)[1],1] <- dim(newcces)[1]
> colnames(table3) <- c("Posterior Mean", "Posterior SD", "Best Nonzero HPD")
> print(xtable(table3, digits = 3), type = "html", file = "table3.html")
> 
> 
> 
> melted_BZ_labeled <- melt(BZ_data)
No id variables; using all as measure variables
> 
> #melted_data[grep("*[6]",melted_data$variable),]$value <- -melted_data[grep("*[6]",melted_data$variable),]$value
> 
> melted_BZ_labeled <- melted_BZ_labeled[(1:nrow(melted_BZ_labeled))[-grep("Constant", melted_BZ_labeled$variable)],]
> 
> # extract results dealing only with personality
> melted_BZ_B5 <- melted_BZ_labeled[grep(paste(personality,collapse="|"), melted_BZ_labeled$variable),]
>   
> 
> 
> ### one plot for gb (opinion slopes), ge (saliency intercepts), and gd (decisiveness)
> bigplot.frame <- rbind(data.frame(value = melted_BZ_B5[grep("beta_b", melted_BZ_B5$variable),], class = 1),
+                        data.frame(value = melted_BZ_B5[grep("beta_e", melted_BZ_B5$variable),], class = 2),
+                        data.frame(value = melted_BZ_B5[grep("beta_d", melted_BZ_B5$variable),], class = 3))
> 
> colnames(bigplot.frame) <- c("variable", "value", "class")
> bigplot.frame$class <- factor(bigplot.frame$class, labels = c("Opinion Slopes", "Saliency Intercepts", "Decisiveness"))
> levels(bigplot.frame$variable) <- list(Openness = c("beta_b[Openness]", "beta_e[Openness]", "beta_d[Openness]"),
+                                        Conscientiousness = c("beta_b[Conscientiousness]", "beta_e[Conscientiousness]", "beta_d[Conscientiousness]"),
+                                        Extraversion = c("beta_b[Extraversion]", "beta_e[Extraversion]", "beta_d[Extraversion]"),
+                                        Agreeableness = c("beta_b[Agreeableness]", "beta_e[Agreeableness]", "beta_d[Agreeableness]"),
+                                        Neuroticism = c("beta_b[Neuroticism]", "beta_e[Neuroticism]", "beta_d[Neuroticism]"))
> 
> # figure 6
> bigplot <- ggplot(bigplot.frame, aes(x=variable,y=value)) + 
+     geom_violin(alpha=0.6) +
+     theme_bw() +
+     facet_grid(~ class) + 
+     #  uncomment this line and comment out the one below for a plot with all variables
+     #  scale_x_discrete(labels=varNames[-1]) +
+     scale_x_discrete(labels=varNames[c(2:6,14)]) +
+     theme(axis.title.y = element_blank(), legend.title=element_blank(), legend.position="none", panel.margin = unit(1, "lines")) +
+     ylab("\nCoefficient Estimate") +
+     geom_hline(yintercept = 0, linetype="dotted") +
+     coord_flip() + 
+     scale_fill_manual(values = "white") +
+     stat_summary(fun.y=median,geom='point') + 
+     stat_summary(fun.data=middle50HPD,geom='errorbar',width = 0.075) +
+     stat_summary(fun.y=middle80HPD,geom='line',size = 0.5)
> 
> cairo_ps(file="bigplot.eps", height = 4, width = 8)
> bigplot
> dev.off()
null device 
          1 
> 
> # how much of the PD for the Neuroticism coefficient in the opinion slope equations lies to the left of zero?
> sum(bigplot.frame[which(bigplot.frame$class == "Opinion Slopes" & bigplot.frame$variable == "Neuroticism"),2] < 0)/
+ length(bigplot.frame[which(bigplot.frame$class == "Opinion Slopes" & bigplot.frame$variable == "Neuroticism"),2] < 0)
[1] 0.98421
> 
> # saliency intercept?
> sum(bigplot.frame[which(bigplot.frame$class == "Saliency Intercepts" & bigplot.frame$variable == "Neuroticism"),2] < 0)/
+ length(bigplot.frame[which(bigplot.frame$class == "Saliency Intercepts" & bigplot.frame$variable == "Neuroticism"),2] < 0)
[1] 0.85468
> 
> # and for Decisiveness?
> sum(bigplot.frame[which(bigplot.frame$class == "Decisiveness" & bigplot.frame$variable == "Neuroticism"),2] < 0)/
+ length(bigplot.frame[which(bigplot.frame$class == "Decisiveness" & bigplot.frame$variable == "Neuroticism"),2] < 0)
[1] 0.71194
> 
> 
> 
> ### appendix plots (PNG used because EPS is way too slow)
> plot.labs <- data.frame(Parameter = c(colnames(DKWYG.nofixed[[1]])), Label = c(rep(varNames,4),"Copartisan",candNames[-c(1:2)],"Cutpoint 2", "Cutpoint 3", "Cutpoint 4", "Cutpoint 5", "Cutpoint 6"))
> 
> # opinion intercept
> CairoPNG(file = "Figures/opinion-intercept-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_a", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+               geom_line(alpha = 0.05)
> dev.off()
null device 
          1 
> 
> 
> CairoPNG(file = "Figures/opinion-intercept-density.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_a", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free")
> dev.off()
null device 
          1 
> 
> 
> CairoPNG(file = "Figures/opinion-intercept-running.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_a", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+             geom_line(alpha = 0.3)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> # opinion slope   
> CairoPNG(file = "Figures/opinion-slope-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)           
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_b", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+               geom_line(alpha = 0.05)
> dev.off()
null device 
          1 
> 
> CairoPNG(file = "Figures/opinion-slope-density.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_b", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free")
> dev.off()
null device 
          1 
> 
> CairoPNG(file = "Figures/opinion-slope-running.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_b", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+             geom_line(alpha = 0.3)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> # decisiveness slope    
> CairoPNG(file = "Figures/decisiveness-slope-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                     
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_d", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+               geom_line(alpha = 0.05)
> dev.off()              
null device 
          1 
> 
> 
> CairoPNG(file = "Figures/decisiveness-slope-density.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_d", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free")
> dev.off()
null device 
          1 
> 
> 
> CairoPNG(file = "Figures/decisiveness-slope-running.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_d", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+             geom_line(alpha = 0.3)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> # saliency intercept    
> CairoPNG(file = "Figures/saliency-intercept-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                               
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_e", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+               geom_line(alpha = 0.05)
> dev.off()
null device 
          1 
> 
> CairoPNG(file = "Figures/saliency-intercept-density.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_e", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free")
> dev.off()
null device 
          1 
> 
> 
> CairoPNG(file = "Figures/saliency-intercept-running.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_e", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             facet_wrap(~Parameter, ncol = 3, scales = "free") + 
+             geom_line(alpha = 0.3)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> # saliency slope     
> CairoPNG(file = "Figures/saliency-slope-traceplot.png", height = 2.5, width = 7.5, units = "in", dpi = 300)                                        
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_s", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               geom_line(alpha = 0.05) + 
+             facet_wrap(~Parameter, scales = "free")
> dev.off()              
null device 
          1 
> 
> CairoPNG(file = "Figures/saliency-slope-density.png", height = 2.5, width = 7.5, units = "in", dpi = 300)
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_s", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue')  + 
+             facet_wrap(~Parameter, scales = "free")
> dev.off()              
null device 
          1 
> 
> 
> CairoPNG(file = "Figures/saliency-slope-running.png", height = 2.5, width = 7.5, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "beta_s", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             geom_line(alpha = 0.3) + 
+             facet_wrap(~Parameter, scales = "free")
> dev.off()
null device 
          1 
> 
> 
> 
>               
> # candidates    
> CairoPNG(file = "Figures/stimuli-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                                                  
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "gamma", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               facet_wrap(~Parameter, ncol = 1, scales = "free") + 
+               geom_line(alpha = 0.05)
> dev.off()              
null device 
          1 
>            
> CairoPNG(file = "Figures/stimuli-density.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "gamma", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue') + 
+             facet_wrap(~Parameter, ncol = 1, scales = "free")
> dev.off()
null device 
          1 
> 
> CairoPNG(file = "Figures/stimuli-running.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "gamma", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             facet_wrap(~Parameter, ncol = 1, scales = "free") + 
+             geom_line(alpha = 0.3)
> dev.off()
null device 
          1 
> 
> 
>               
> # cutpoints  
> CairoPNG(file = "Figures/cutpoints-traceplot.png", height = 10, width = 10, units = "in", dpi = 300)                                                              
> ggs_traceplot(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "tau", burnin = F)) + 
+               theme_bw() + 
+               ylab('Value\n') + 
+               xlab('\nIteration') + 
+               facet_wrap(~Parameter, ncol = 1, scales = "free") + 
+               geom_line(alpha = 0.3)
> dev.off()              
null device 
          1 
> 
> CairoPNG(file = "Figures/cutpoints-density.png", height = 10, width = 10, units = "in", dpi = 300)              
> ggs_density(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "tau", burnin = F)) + 
+             theme_bw() + 
+             ylab('Density\n') + 
+             xlab('\nValue') + 
+             facet_wrap(~Parameter, ncol = 1, scales = "free")
> dev.off()            
null device 
          1 
> 
> CairoPNG(file = "Figures/cutpoints-running.png", height = 10, width = 10, units = "in", dpi = 300)
> ggs_running(ggs(DKWYG.nofixed, par_labels = plot.labs, family = "tau", burnin = F)) + 
+             theme_bw() + 
+             ylab('Value\n') + 
+             xlab('\nIteration') + 
+             facet_wrap(~Parameter, ncol = 1, scales = "free") +
+             geom_line(alpha = 0.3)
> dev.off()
null device 
          1 
